% estimate impulse responses to government spending shock 
% recursive identification
% this version: May 2011

clear all; clc; close all;

load us_quarterly       % generated by getdata.m

saveirf      = 1;

rand('state',1);
nptsvar     = 30;       % number of points of IRF in VAR
teststatio  = 0;        % tests wether VAR is explosive
confidence  = 1.645;    % (1.645: 90 percent); 
ndraws      = 1000;     % number of replications in bootstrap
nlagq       = 4;        % number of lags
beg_samp    = (82-47)*4+1; 
end_samp    = (107-47)*4+4;

disp(['Sample: ' num2str(sampleq(beg_samp+nlagq)) ' - '  num2str(sampleq(end_samp)) ' (dependent variables)  '  ])
gy = mean(g(beg_samp:end_samp)./y(beg_samp:end_samp));
cy = mean(c(beg_samp:end_samp)./y(beg_samp:end_samp));
iy = mean(inv(beg_samp:end_samp)./y(beg_samp:end_samp));


disp(['importshare: ' num2str(mean(imp(beg_samp:end_samp)./y(beg_samp:end_samp)))])
disp(['spending-share: ' num2str(mean(g(beg_samp:end_samp)./y(beg_samp:end_samp)))])


for jrun = 1:3;

    if jrun==1;
        namshort = char('g','y','c','lr','rx','dy','inf');       
        xdata = [log(g) log(y) log(c) lrr log(rx) dy  inflation];   
        yshare = [gy 1 cy 1 1 1 1];
    elseif jrun ==2;
        namshort = char('g','y','i','lr','rx','dy','inf');       
        xdata = [log(g) log(y) log(inv)  lrr log(rx) dy inflation];
        yshare = [gy 1 iy 1 1 1 1];
    elseif jrun ==3;
        namshort = char('g','y','nx','lr','rx','dy','inf');       
        xdata = [log(g) log(y) nx lrr log(rx) dy inflation];
        yshare = [gy 1 1 1 1 1 1];
    end

    %% specify deterministics
    [obs nvar]=size(xdata);
    exo = []; 
    for ik=1:obs; ltrend(ik)=ik; end
    for ik=1:obs; exosq(ik)=ik^2; end
    exo = [ones(obs,1) ltrend'];
    exo = exo(beg_samp:end_samp,:);
    samplelength=length(xdata(beg_samp:end_samp,:));

    %% estimate impulse responses
    [impzout,erzout,azeroout,a0betazout,var_explode]=...
        estimateIRF(xdata(beg_samp:end_samp,:),nlagq,nptsvar,teststatio,exo);       
    irf1 =  impzout(:,:,1);
    
    %% bootstrap IRF
    simdout=simdata(xdata(beg_samp:end_samp,:),...
    nlagq,nptsvar,ndraws,impzout,erzout,azeroout,a0betazout,exo);
    simimpzout=[];
    for zx=1:ndraws;
        [simimpzout,erzout,azeroout,a0betazout,var_explode] = estimateIRF(simdout(:,:,zx),...
        nlagq,nptsvar,teststatio,exo);                   
        irf1_mc(zx,:,:) = squeeze(simimpzout(:,:,1));
    end
    irfse1 = squeeze(std(irf1_mc(:,:,:),0,1));
    

    %% plot
    figure
    time = [0:nptsvar-1];
    for zx=1:nvar;
        subplot(3,3,zx)
        mycolors = [1 1 1; .85 .85 .85];            
        colormap(mycolors);

        rescale = max(irf1(1:10,1))*(yshare(1)/yshare(zx));
        point = irf1(:,zx)./rescale;
        sepoint = irfse1(:,zx)./rescale;
        ub = point+confidence*sepoint;
        lb = point-confidence*sepoint;  

        hold on;   
        gx=[lb ub-lb];
        area(time,gx,min(lb),'EdgeColor','none');
        plot(time,zeros(nptsvar,1),'k--','LineWidth',1);
        plot(time,point,'b-','LineWidth',1);
        axis([-.2 nptsvar-.8  min([1.5*min(lb) 0]) 1.5*max(ub)]);box on;       
        box on;
        title(deblank(namshort(zx,:)))        
    end
    
    if saveirf
        
        disp('writing figures to files...')
        for zx=1:nvar;            
            figure
            set(gcf,'Visible','off')
            mycolors = [1 1 1; .85 .85 .85];            
            colormap(mycolors);

            rescale = max(irf1(1:10,1))*(yshare(1)/yshare(zx));
            point = irf1(:,zx)./rescale;
            sepoint = irfse1(:,zx)./rescale;
            ub = point+confidence*sepoint;
            lb = point-confidence*sepoint;  

            hold on;   
            gx=[lb ub-lb];
            area(time,gx,min(lb),'EdgeColor','none');
            plot(time,zeros(nptsvar,1),'k--','LineWidth',1);
            plot(time,point,'b-','LineWidth',4);
            axis([-.2 nptsvar-.8  min([1.5*min(lb) 0]) 1.5*max(ub)]);box on;       
            set(gca,'FontSize',40)
            box on;
            saveas(gcf,['..\figures\' num2str(jrun) '_' deblank(namshort(zx,:)) '.eps'],'psc2')          
        end

    end
end